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{—} | Abstract Recently observed oscillations in the solar atmosphere have been inter- 

£*\j ■ preted and modeled as magnctohydrodynamic wave modes. This has allowed the 

estimation of parameters that arc otherwise hard to derive, such as the coronal 
magnetic-field strength. This work crucially relies on the initial detection of the 
oscillations, which is commonly done manually. The volume of Solar Dynamics 

>0 . Observatory (SDO) data will make manual detection inefficient for detecting all 

of the oscillating regions. An algorithm is presented which automates the detec- 
tion of areas of the solar atmosphere that support spatially extended oscillations. 

jvS | The algorithm identifies areas in the solar atmosphere whose oscillation content is 

described by a single, dominant oscillation within a user-defined frequency range. 
The method is based on Bayesian spectral analysis of time-scries and image 
filtering. A Bayesian approach sidesteps the need for an a-priori noise estimate to 
calculate rejection criteria for the observed signal, and it also provides estimates 

-^ . of oscillation frequency, amplitude and noise, and the error in all these quantities, 

jrt ' in a self-consistent way. The algorithm also introduces the notion of quality mea- 

sures to those regions for which a positive detection is claimed, allowing simple 
post-detection discrimination by the user. The algorithm is demonstrated on 
►> . two Transition Region and Coronal Explorer (TRACE) datasets, and comments 

l/"} ' regarding its suitability for oscillation detection in SDO arc made. 

t^. ! 

qv , Keywords: Sun: active region, Sun: magnetic field 

o 



1. Introduction 

The Solar and Relio spheric Observatory (SOHO) and TRACE missions es- 
tablished that the solar corona supports observable oscillations. The extensive 
literature on the theory of oscillations in coronal flux tubes has been used 
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in conjunction with observationally derived parameters to deduce conditions 
in the corona. The study of these oscillations is called coronal seismology (in 
analogy to terrestrial seismology). Their general features, as determined by 
current analyses, and their implications for the physics of the corona are sum- 
marized briefly below; more detailed reviews can be found in lDe Moortell (I2005P 
andTNakariak ov and Verwic hte (2006). Quasi-periodic oscillations, interpreted as 
propagating slow waves have been detected in coronal plumes SOHO/Ultraviolct 
Coronal Spectrograph fUVCS: lOfman et al\\\M7) : SOHO/Extremc Ultraviolet 
Imaging Telescope (EIT: Deforest and Gurmanl I1998|) . A similar phenomenon 
(at higher frequency) has also been detected at the base of coronal loops (Berghmans and Clette|| 



ITgg§l|Nightingale, Aschwanden, and Hurlburtlrrg<M|De Moortel, Ireland, and Walsh[ | 
120001 TRobbrecht et qZ.l I2001|) . Thev are interpreted as propagating slow waves 
since they travel at approximately the sound speed and are seen as intensity 
(and therefore density) variations propagating and decaying away from the base 
of the flux tube through the corona. The observed periodicities fall into distinct 
ranges around three and five minutes (jDe Moortel et q771 12002)) . suggesting that 
the oscillations are due to different connectivity to either sunspot or transition- 
region moss magnetic field respectively. Coronal seismological applications for 
these oscillations include determination of the connectivity of the photosphere 



to the corona (de Pontieu, Erdelyi, and De Moortel I2005|) as well as deriving 



information on the coronal heating function since the wave observables (e.g., 
period and damping length) are strongly dependent on the thermal conditions 
of the corona (|De Moortel and Hood! I20U51 IDe Moortel and Hood! I2UM)) . 

All of these studies have in common that the oscillating feature was discovered 
by a manual examination of data in regions known (or suspected) to contain 
oscillating material. Although a successful strategy for detecting oscillations, it 
is clearly impractical for ever-increasing rates of data acquisition. For example, 
the SDO mission is acquiring £l.4 Tb of science data per day, around 1000 
times the data rate of SOHO. Much of this data will be taken at cadences 
and spatial resolutions comparable to the best that TRACE can produce, the 
key difference from TRACE being that SDO will provide a near continuous, 
full disk coverage of the Sun in multiple wavelengths simultaneously. Given the 
diagnostic possibilities of these waves, and the large amount of data that will 
become available, an automated detection algorithm is necessary. There are at 
least five published algorithms, which are briefly described below. 



Nakariakov and K ing (2007) use a thresholded fast Fourier Transform to find 
locations in TRACE data that may support an oscillatory signal. The threshold 
level is defined as three -four times the average FFT power; if the maximum 
FFT power is above this level then the frequency at which that power occurs 
is assumed to be real. Since this method relies on the FFT, it is very fast. 
The definition of the threshold is chosen to speed up the algorithm in compar- 
ison to calculating a threshold; for example, implementing the randomization 
method of ILinnell-Nemec and Nemecl (|1985p requires at least an estimated 150 
times more FFT calculations. ILinnell-N emec a nd Nemecl (|1985p establish fre- 
quency acceptance/rejection criteria based on probabilistic arguments whereas 



Nakariakov and King (200TJ) use essentially empirically derived arguments to 



establish the threshold. Further, identification of contiguous groups of pixels 
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that might form an oscillatory region is carried out by manual inspection. To 
aid identification, the authors assume a null hypothesis that the entire time- 
series consists solely of Gaussian distributed noise, and so the chance that any 
two neigboring pixels contain the same frequency is small. Hence, if two or more 
neighboring pixels do satsisfy the selection criteria, then it is likely they are 
physically connected. These assumptions ignore two effects: firstly, neighboring 
pixels are not statistically independent due to the point spread function of the 
instrument. Secondly, the Sun is observed to have physically connected struc- 
tures, such as loops, that exist over many pixels in at least one direction, and 
so by continuity, one can expect that neighboring pixels do influence each other. 
Hence the assumption that neighboring pixels are independent of each other 
is an over-simplification of the nature of the image. The authors suggest that 
the algorithm be used to identify regions in the data worthy of further study, 
although there is no quoted method of automated region identification. 

IDe Moortel and McAteerl ([2004) describe an automated oscillation detection 



algorithm based on the wavelet analysis routines of Torrence and Compo (I1998[) 
and the analysis procedure of llreland et al.\ (|1999|) . The algorithm finds signif- 
icant wave packets ranging from single to multiple wave cycles in duration, by 
a wavelet power/confidence level comparison against the null hypothesis that a 
given time-series is Gaussian distributed noise. Shorter duration detections are 
rejected. Contiguous regions of multi-cycle duration wave packets are found in 
the data, but are identified and isolated manually by inspection. 

Syc h and Nakariakov|(|200"S|) base their detection algorithm on pixclized wavelet| 
filtering (PWF) of a three-dimensional data cube (x,y,t). This too is based on 



the wavelet analysis routines of Torrence and Compo| (|1998[) . but has a more 



complex treatment of the resultant wavelet spectrum. Regions of interest are 
found by first calculating a variance map (jGrechnevi I2003P of the signal; regions 
with high variance are candidate oscillatory regions (note that this must also 
imply the removal of a background trend in order for the variance to measure 
an oscillatory signal, and not the trend). Wavelet spectra are calculated for 
those pixels, and only the "significant" pixels are retained (what constitutes a 
significant signal is not stated explicitly). The routine analyzes further (than 
Do Moortel and McAtccr, 2004) the temporal evolution of the oscillation and so 
can differentiate between standing and traveling waves. 

The algorithm presented by |McIntosh, de Pontieu, and Tomczyk] (|2008p can 
also differentiate between standing and traveling waves. The algorithm begins 
by Fourier transforming the entire data cube, performing cross correlations 
with neighboring pixels in narrow frequency bands, and filtering the results (by 
thresholding on various quantities, such as eliminating areas where the relative 
error in the calculated phase speed is large) to determine groups of pixels that 
are highly correlated in both time and space. This correlation technique allows 
the discrimination of standing and traveling waves, and the calculation of other 
parameters such as the phase speed and the propagation angle. There is an 
implicit null hypothesis in the algorithm: at one stage, only time-series with a 
high coherence are accepted for further analysis. The null hypothesis here is that 
the candidate time-series is pure noise, with the additional assumption that pairs 
of noisy time-series have low coherence, and so can be rejected. However, this is 
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not strictly true; IChatfieldl (|1996p shows an example of two different noisy time- 
series that have a perfect coherence over all spectral frequencies because both 
time-series are generated from the same noise process. Highly coherent patches 
of non-oscillatory material are probably filtered out in the next stage of the 
algorithm, where regions with a poorly determined phase speed (large error) are 
discarded; this, however, has not been explicitly tested (S. W. Mcintosh, private 
communication, 2008). 

It is clear that there is an increasing amount of effort in finding oscillatory 
regions and identifying waves in the solar atmosphere. All the above algorithms 
show promising results and avenues for further work. The algorithm presented 
here seeks to find regions in the data which support oscillations via Bayesian 
time-series analysis. This involves calculating explicitly the probability that a 
time series supports an oscillation of a given frequency. This is in distinction to 
the methods above, which rely on statements about null hypotheses in order to 
determine if an oscillation is present. Section[5]introduces Bayes' Theorem and an 
application of it to time-series analysis. Section[3]describes a detection algorithm 
based on the results of Section [21 whilst Section @] describes the application of 
this algorithm to some example datasets from the TRACE mission. Finally, 
Section |B] discussions some further applications of Bayesian time-series analysis 
and automated EUV detection algorithms. 



2. Bayesian Time-Series Analysis 

Denoting by p(a\b) the conditional probability that proposition a is true, given 
that proposition b is true, Bayes' theorem is 

mDtI) = mD*mii (1) 

p(D\I) 

where H is the hypothesis to be tested, D is the observation, and / is any 
applicable prior information that we have before making the observation ( |Bayes| 
I1763p . The left hand side p(H\D,I) is called the posterior probability of the 
hypothesis, given the data and the prior information, and it encapsulates the 
available knowledge about the hypothesis. The quantity p(H\I) is called the prior 
distribution and represents what we know about H prior to any data collection. 
Often a prior describes a probability distribution of likely parameter values. The 
sampling distribution or likelihood (p(D\H,I)) represents the likelihood of the 
data given the hypothesis (as well as any prior information). The quantity p{D\I) 
is the prior probability of the data; it is absorbed into a normalization constant, 
and does not affect the following analysis for a given model. Equation ([l) is 
very general and is not restricted to the mathematical equations: any logical 
proposition can be treated in a Bayesian context ( |Jaynes[ 120031 |Gregoryj I2005P . 

2.1. Signal Containing a Single Frequency 

For the purposes of detecting oscillations in solar time-series d(ti), the linear 
model 

d{ti) = b\ cosLoti + &2 sin(wti) + Xi (2) 
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is used. This assumes that all time-series are modeled as a single oscillation plus 
Gaussian distributed noise. Under some simplifying assumptions (N >> 1 and 
that there are no low-frequency oscillati ons present), |Jaynes| (|1987p . IBretthorstl 



(|1988p . and O Ruanaidh and Fitzgerald ( 1996) derive an expression for the prob 



ability that the time-series contains an oscillation of frequency uo for the model 
oscillation above, that is, 



p{u\d,I) oc 



where p = l/^X^i ^? an d 



C(w) 



N 



2CV) 
up 1 



E* 



(2-n)/2 



(3) 



(4) 



is the Schuster periodogram (|Schustertll898[) . The full details of the derivation of 
Equation ([3]) are contained in the appendix. This is a probability measure that 
a given frequency is present in the data that does not require explicit knowledge 
of the Gaussian noise present. Equation ([3]) can analyze unevenly sampled data 
due to its use of the Schuster periodogram. If the analysis frequencies are 



j p = 2irp/n, < p < n/2 — 1, 



(5) 



Equation ^ is proportional to the power of the fast Fourier transform (FFT) 
of the data ( |Cooley and Tukey[ ITM51 IGhatfieldl IT9961) . The speed of the FFT 
transform can be exploited for data which are evenly or close to evenly spaced. 



2.2. Example Analysis 

Figure Q] demonstrates the application of the above equations to the analysis of 
a single time-series. The true oscillation is an evenly sampled simple cosine at 
frequency / = 3.33 mHz (100 data points at a cadence of 26.4 seconds) with unit 
amplitude. The fake data is this oscillation corrupted by Gaussian distributed 
noise with deviation of a = 1, a signal-to- noise ratio of 1 (Figure |T^a)). The 
probability function, Equation [31 is highly peaked close to the true frequency 
and the probability that the signal contains a frequency within the range 3.0- 
3.5 mHz exceeds 99% (Figure [TJb)). Figures QJc,d) snow the measured basis 
function amplitudes and Gaussian deviation estimate respectively as a function 
of frequency. Both lie close to their true values at the true frequency, as indicated 
by the solid vertical line. The quoted values shown in each plot (Figurc[ljb,c,d)) 
are found at u = ui, the most probable frequency. An error estimate is given by 
the standard deviation 



Nu 



sd 



= ^2p(wi)[<f>(u 



mr 



(6) 
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Figure 1. An example analysis using the Bayesian approach of Section [2] Panel (a) shows the 
true oscillation, the observed (noisy) data, and the estimate. Panel (b) shows the probability- 
density function arising from the observed data. Panel (c) shows the estimated amplitudes as 
a function of frequency [/] , and panel (d) shows the estimated noise as a function of frequency. 
Final quantity estimates and an estimated error are shown also. See Section [2] for more detail. 



where (j){uj) stands for the analyzed frequency range, the distribution of ampli- 
tudes (found from Equation (fT5j) - see appendix) and the distribution of a (found 
from Equation (|15p - see appendix) as a function of the analysis frequency set 
Wj, 1 < i < N^ (such as the Fourier set, Equation ([5])). 

These estimates of frequency probability, oscillation amplitude and Gaussian 
noise deviation have all been obtained without explicitly performing a least- 
squares fit, or with no special a-priori knowledge of the noise in the signal other 
than the assumption of a Gaussian distribution. 

The peakcdness of the Bayesian probability-density function immediately 
suggests that a search for portions of the {to} parameter space that contain 
most of the probability are the values that are of greatest interest. In addition, 
prior knowledge from previous studies of solar atmospheric oscillations suggests 
which regions of the parameter space are of interest. These two observations are 
combined below to take the first steps towards creating an automated oscilla- 
tion detection algorithm which can identify oscillation regions and return useful 
information, such as their amplitude. 
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3. Detection using Bayesian Spectral Analysis 

3.1. Data 

The data that we will use are three dimensional datacubes (x,y,t), with the x 
and y locations referring to spatial locations (pixels) on the Sun, and t, time, 
time-series are formed by choosing a particular location on the Sun and extract- 
ing a one-dimensional time-series (occasionally super-pixels formed by the sum 
of neighbouring pixels are used to increase the signal-to- noise ratio). 

3.2. Probability and Frequency Bands 

In order to use Equation ([3]) to detect oscillations, further assumptions must be 
made. The algorithm assumes that the range of possible frequencies is limited to 
that spanned by the Fast Fourier Transform applied to a time-series of similar 
length. This permits us to normalize the probablity density function over a fixed 
range of frequencies. We are commonly interested in detecting oscillations in 
given frequency bands (say the three or five minute oscillation frequency bands), 
and so we calculate the probability that the oscillation in a given pixel lies within 
the range [wi,^], that is, 

Pwi,w a = / p(u)\D,I)du) (7) 

at every point in the image. Large values of p Ul M2 indicate that the true oscil- 
lation frequency is very likely to be within the range [wi,^]- 

The Bayesian formulation handily yields both computational and logical ad- 
vantage over other frequcntist approaches for oscillating-pixel detection. The 
prime derived data product for each pixel is a probability distribution describ- 
ing the probability of a given frequency in the data and so there is no need 
to perform secondary, often computationally expensive calculations, such as 
randomization tests (jLinn cll-Nemcc and Ncmcc, 1985; O'Sh ca et ail I2001[) to 
assess the probability that the frequency is present or not. Frcquentist-bascd 
approaches to detecting oscillations rely on calculating the expectated value 
of the Fourier transform power at a given confidence level, assuming that the 
observed data arises from the null hypothesis that the time-series is pure noise. If 
the Fourier transform power of the observed time-series exceeds the expectation 
value, then the null hypothesis is rejected at that confidence level. This, however, 
strictly cannot be used to imply the presence of an oscillation; we have merely 
rejected the null hypothesis. In comparison, Equation ([7]) directly calculates the 
probability that the frequency of oscillation lies in the range [wi , W2] , under the 
modeling assumptions of Section [5J 

3.3. An Algorithm 

Figure [5] describes the general-purpose algorithm used to find oscillating spa- 
tial locations in the data cube. Step 2 relies on instrument calibration routines 
provided by instrument teams, and on standard solar de-rotation routines, both 
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Algorithm: Bayesian probability based oscillation detection 

1 begin 

2 Prepare data: apply instrument calibrations and de-rotate data cube. 

3 Dctrend data and ensure data has zero mean. 

4 for each time-series d x ,y{t) in the data cube (x, y, t) 

5 Calculate p Ul M2 for given Wi,W2 a t [x,y). 

6 end 

7 Generate map M of the spatial distribution of probability p Ul )U2 

8 Filter map to find the "highly probable" oscillation areas. 

9 Report these areas and derive useful parameters. 

10 end 



Figure 2. Pseudo-code algorithm to find oscillating locations in (x,y, i) data cubes. 

provided in the IDL/Solarsoft package. Time-series de-trending (Step 3) is nec- 
essary for two reasons. The first is dictated by our choice of oscillation model, 
Equation @. This model assumes that the time series contains a single oscilla- 
tion only, and no other features. (It is certainly possible to define other models 
that do contain a background trend, or multiple frequencies, and calculate prob- 
ability density functions for those time-series. However, these more sophisticated 
analyses are not necessary for us to make progress in the current application of 
locating oscillating material). Secondly, strong background trends can pollute the 
Fourier power spectrum with spectral power unrelated to the oscillation. This 
can lead to the mis-identification of peaks in the power spectrum as oscillatory 
when they are not. It should be noted that the influence of background trend on 
locating oscillations of a given frequency will influence all proposed automation 
algorithms. 

De-trending in Step 3 is accomplished by subtracting a running average of 
the time-series (a window of size R seconds is slid across the time-series and 
the running average of the data lying entirely within that window is calculated). 
Since previous experience has told us that there are oscillations of interest that 
have periods less than 500 seconds, it seems to us that 500 seconds is a reasonable 
to chose for these data (jAsc hwandc n et ai[ [2002) . Timescales longer than R can 
be considered to be associated with the background trend. Timescales shorter 
than R potentially support an oscillation. Fourier power spectra of smooth time- 
series that have had their background trend subtracted have extremely low power 
at frequencies less than that corresponding to R. Therefore, these low frequencies 
corresponding to the background trend have an extremely low probability and 
are not selected. 

Step 5 finds the probability that the time-series has a single frequency in the 
range [ui, uj2\. There are several ways of defining [ui, U2\. In the work below the 
algorithm begins by first finding the location Wmax of the highest peak in the 
probability distribution function (PDF) p(u>\D,I). If [fii < Wmax < ^2], where 
[f2i, fy arc defined by the user (a frequency filter), then the algorithm proceeds 
by stepping away from the peak to find the nearest turning points in the PDF, 
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Figure 3. Example image data and probability and frequency maps. Panels a, b, and c refer 
to TRACE 171A data taken on 1st July 1998 and analyzed for oscillations in the three-minute 
frequency band 4.17 < /max < 8.33 mHz), and panels d, e and f refer to TRACE 171A 
data taken on 14 July 1998 and analyzed for oscillations in the five-minute frequency band 
(2.78 < /max < 4.17 mHz) . 



located at [wi,^]- The probability [p Ul ,ui 2 ] IS then calculated. This creates the 
probability map M of Step 7. 

Figure |3Ia,b,c) shows some example data, and its analysis up to Step 7. Two 
datasets are examined, both in the TRACE 17lA wave-band. The first dataset 
consists of 391 observations with a mean 31s cadence taken on 1 July 1998 
12:03:10 UT - 15:28:52 UT. The resulting probability map and its concomitant 
frequency map are shown in Figures [3^b,c) respectively. There appears to be 
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a single large area that appears to contain a significant oscillation at differing 
frequencies and probabilities. 

Cleanly extracting oscillating features is much more difficult for the second 
example dataset shows in Figures [3][d,e,f). Here there are many groups of pixels 
having similar frequencies and high probability close to each other. In addition, 
there arc many pixels that have high probability scattered all over the field of 
view, further complicating the task of extracting a single group. The eye is very 
good at picking out such groups, and the algorithm attempts to mimic some 
of this behavior. The detection algorithm attempts to mimic this process by 
examining the local spatial distribution of probability. Firstly, it is assumed that 
nearby pixels of high probability are physically related, even if their frequencies 
may be different. However, these groupings may not be contiguous due to noise, 
for example. This may be overcome by smoothing the probability map at differ- 
ent lcngthscalcs, and flagging those areas which exceed a certain threshold. This 
has the effect of the allowing discontinuous areas of high probability to merge, 
in a similar way to how the eye may integrate physically "close" oscillating areas 
into one group. More formally, Step 8 implements 

1. for m an integer in the range [Xi,!^] : 

a) Smooth the probability map [M] with a Gaussian filter of width m. 

b) Generate a mask [Z m ] locating all the areas in the smoothed map that 
have a smoothed probability at or over 0.5 

The lengthscales m allow the user to examine the average spatial probability 
structure of the data on multiple lengthscales. This procedure is meant to 
mimic the way the eye examines such probability maps, where the eye finds 
it easy to integrate over neighbouring high probability features. The upper 
lengthscale [L2] describes the maximum distance that two pixels are assumed 
to be potentially physically related to each other. Smaller maximum values 
of m cause the smoothed probability masks to break up more, whereas larger 
values can cause the creation of groups of pixels which are essentially unre- 
lated. The masks [Z m ] found in step lb denote the areas that are more likely 
than not to support contiguous oscillations over the lengthscale m. The net 
effect of step 1 is to implement a simple multi-scale analysis of the spatial 
probability structure. 

2. Add all the masks [Z m ] together. This is the mask — Z — of all the candidate 
pixel groups that may support a significant oscillations. 

3. Remove from Z all the candidate pixel groups that form a group smaller than 
g pixels. Nominally, we set a minimum group size g = L\\ that is, we look 
for groups of pixels that are connected over the largest distance set in the 
multi-scale analysis of step 1. 

4. For a candidate group of area A, remove that group from Z if it has more 
than hA zero-probability pixels, < h < 1. The remaining pixel groups are 
considered to be regions in the solar atmosphere that support oscillations in 
the range f2i, f?2- 

In the results below, L\ = 1 and L2 = 4, g = 16; these are chosen in order 
to allow smaller width strands to be found in the data. The last part of Step 
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Table 1. Quantities reported for a contiguous group [G] of pixels found by the detectionB 

algorithm of Figure [2] 

quantity definition units 

/ average frequency in a group G mHz 

o$ standard deviation of frequency mHz 

A area of G with a non-zero probability area (px) 

Ao.g5 area of G with p wl]W2 > 0.95 area (px) 

F yl/(area of group G) dimensionless 
p average probability of all pixels having a non-zero probability dimensionless 

Q Fp (quality of the group G) dimensionless 

E AQ (the equivalent area of the group G) area (px) 



8 acknowledges that in performing the probability-map smoothing, pixels that 
contain no oscillation in the range Oi , O2 can be swept into a candidate pixel 
group A, and this must have a limit lest the candidate pixel group have too many 
lacunae. The final step (Step 9, Figurc[2]) is to report properties of the remaining 
pixel groups such as the average frequency [/] and the standard deviation [ay] of 
the frequency. The full table of reported results is given in Table [TJ In any given 
analysis, the regions found result from a series of filters and thresholds made by 
the user. However, not all of the regions that survive this process are necessarily 
equivalent, they have merely satisfied some set of criteria. These quantities 
attempt to measure differences between the surviving oscillatory regions, and 
allow the user to further discriminate them at their discretion. Clearly the size 
[area A] of the pixel group is important, and the number of highly probable pixels 
[A0.95]. The quantity F measures how complete the group is, and p measures the 
average probability that the average oscillation frequency of G is indeed between 
Sli , r^2- The quality [Q] is an average probability for the entire group, and the 
equivalent area [E] is the area that the group would have if it were entirely 
complete (no lacunae) and every oscillation detected in it had p u>1 . u>2 = 1. These 
measures, along with maps of the detected regions showing frequency, amplitude, 
and noise estimates (along with error estimates to each of these quantities), are 
the final derived products of this automated detection analysis. 



4. Results 

Data was analyzed from two TRACE observing sequences. Data was prepared 
using standard IDL/Solarsoft TRACE routines (TRACE.PREP with the key- 
words /wave2point, /unspike, /destreak, /deripple, /norm, /float switched on) 
and derotated using IDL/Solarsoft SHIFT_XY to calculate the motion of the 
observed piece of Sun, and cubic two-dimensional interpolation to shift subse- 
quent images back to the location of the first one. In addition, a strip of data 
of width ten pixels at the extreme right of the image is removed, as this is 
the location of severe image distortion due to image derotation. Further, all 
time-series are detrended by removing a running average taken over aii= 500 
second timescale (see Step 3, Figure [5]). In the following results, signal-to-noisc 
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.2 ~ 2 

ratio (SNR) is calculated as the estimated amplitude y bi +62 (see Equation 
©) divided by the Gaussian noise estimate a. 

4.1. 1 July 1998 

The original data was taken on 1 July 1998 12:30:01 - 14:24:41 UT with an 
average cadence of 31 seconds (220 samples) and an image size of 512 x 512 
pixels of equivalent size 0.5". Data arc also 2x2 summed in space to increase 
signal-to-noise ratio (this has the consequent effect of reducing the algorithm 
run time). The results of searching in a wideband three-minute range (120 - 240 
seconds, or 4.17-6.33 mHz) are shown in Figures @] (171 A) and[5](195 A) and 
in Table [U(a,b). Most of the field of view does not contain material oscillating 
in the analysis frequency range. Only one group of pixels survives the filtering 
process, at the base of a coronal- loop fan (Figure E]). However, it is interesting 
to note that there are areas of the probability maps Figures HJb),[5Jb) that have 
effectively zero probability of supporting an oscillation in the frequency band, 
whilst others have a non-zero and low-probability We speculate that with the 
improved SNR of SDO, many of these low probability oscillations will become 
much more likely, leading to the detection of many more spatially distributed 
signals in the data. 

The single oscillating pixel group found (Figure [5]) is in the same location as 
that identified manually by King et a l. (2003) and automatically by Mcintosh, de Pontieu, and Tomczykj 



(2008J. The region is distinctly different from others in the data. Figures 0Jf) 
and[5]Jf) show that the detected oscillations in this pixel group has a fairly low 
estimated signal-to-noise ratio. The estimated error in the frequency decreases 
with increasing SNR, as expected. 

The results of searching in the five- minute frequency band (240-360 seconds, 
or 2.78-4.17 mHz) arc shown in Figures [7] (171 A) and [5] (195 A) and in 
Table [5Jc,d). The probability maps Figures HJb), EJb) look very different in 
this frequency range, with low-probability oscillations scattered across the field 
of view. It is noticeable, however, that the higher-probability oscillations are 
concentrated in the core of the active region, in the regions that appeared empty 
in Figures HJb) , EK^) ■ The algorithm does qualify some pixel groups as support- 
ing oscillations. |King et~aL (|2003|) do not claim any detections in these areas 



(although it is not clear if they looked), and |Mc!ntosh, de Pontieu, and Tomczyk| 
(2008]) find only one coherent pixel group in the same general region as those 
found here. Some of the regions do overlap, suggesting co-temporal and co- 
spatial propagation of oscillations in two different layers of the atmosphere. In 
the image, the areas in question resemble TRACE moss ( |Berger et al.\ [1999; 
Ide Pontieu et ail 119991 IFletcher and de Pontieu! I1999p . However the presence 
of groups of oscillations at multiple temperatures argues that the algorithm has 
found examples of leakage of five minute oscillations from lower down in the 
atmosphere to the upper layers, as described by de Pontieu, Erdelyi, and James| 
(|2004p and |de Pontieu, Erdelyi, and De Moorte | (|2005[) and references therein. 
Finally, the plotted points in Figures[7Jf),[HKf) appear to lie on horizontal lines 
across the plot, which are simply the analysis frequencies. The foregoing analysis 
can be done with many more frequencies (sec Section ^. 1|) , at the expense of using 
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Table 2. Quantities reported (see Tablef^for the automatically detected regions found 
in 1 July 1998 TRACE 17lA and 195A data - see Figures [3] [5] and E] for maps 
and plots of the detected regions listed in (a), (b), (c) and (d) respectively above. The 
quantities Q and E are also listed with their rank [r] when compared to all other regions 
found in the same dataset. 



data 



/±»/ 



Ao.95 



Q [r(Q)] E [r(E)} 



3 min 
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25 
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0.773 
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(b) 195A 


6.00 ± 0.56 
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1.000 


0.788 
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5 min 


















(c) 17lA 
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0.897 


[3] 


13.457 [10] 
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20 


18 


0.870 


0.965 


0.840 


[8] 


16.791 [7] 
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['] 
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67 


39 


0.817 
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7 
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(d) 195A 
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31 


15 


0.912 
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[6] 
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75 


12 
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0.748 
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56.099 [1] 
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3.166 ± 0.254 


16 


6 


1.000 


0.840 
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[3] 


13.433 [8] 


5 


3.240 ± 0.286 


11 


25 
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0.895 


0.816 


[1] 
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14.606 [7] 
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34 


21 
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[7] 


26.452 [3] 



the FFT to perform the analysis. This would slow performance, but would lead 
to a more precise knowledge of the frequencies present. Since we arc primarily 
interested in detecting the presence of frequencies in a wide frequency band, the 
precision of each frequency detected is not as important as their detection in as 
little time as possible. 

4.2. 14 July 1998 

The original data was taken on 14 July 1998 12:45:19 - 13:42:44 UT with an 
average cadence of 73 seconds (sixty samples) and an image size of 512 x 512 
pixels of equivalent size 0.5". Data are also 2x2 summed in space to increase 
signal-to-noise ratio. TRACE was observing NOAA AR 8270 when a GOES class 
M4.6 flare occurred at around 12:55 UT. 
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Figure 4. Example image data and results from the analysis algorithm for TRACE 171A 
data taken on 1 July 1998 (panel (a)). Panel (b) is the probability that a given pixel supports 
an oscillation within the frequency band indicated in panel (c), in this case, a three- minute 
oscillation frequency band. Panel (c) shows the frquency supported at those pixels. Also in- 
dicated in panels (a-c) are the detected oscillation regions. Note that the region agrees with 
that found manually by |King et qZ. 1 1200311 on the following day. The oscillation region is found 
at the base of a coronal loop structure, as has been found by many authors (De Moortel et al. 
2002). Panel (d) shows a map of estimated amplitude for the detected regions, whilst panel 
(e) shows a map of the signal to noise ratio. Finally, panel (f) shows a panel of the detected 
frequency as a function of the signal to noise ratio. The error bars on the abscissa values are 
shown at one-tenth of their actual size in order to better show the lower error values at higher 
signal-to-noise ratio. Table [2ja) shows the values for the detected region for the parameters 
listed in Table Q] 
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Figure 5. Example image data and results from the analysis algorithm for TRACE 195A 
data taken on 1 July 1998 (panel (a)). Panel (b) is the probability that a given pixel supports 
an oscillation within the frequency band indicated in panel (c), in this case, a five- minute 
oscillation frequency band. Note that the region agrees with that found manually by |King et al.\ 
(2003) on the following day. Panels (d-f) show similar maps to those shown and described 
in Figure [4] Table 0b) shows the values for the detected region for the parameters listed in 
Table[l] 



A movie of the analyzed data cube shows that the first eight samples show 
no sign of any flaring. After that, the flare continues for about another 15- 
18 frames. Finally, the flare dissipates and post-flare loops are observable (20 
frames) . The flare event transfers momentum to the surrounding medium, caus- 
ing loops to oscillate. In addition, the whole region evacuates - the flare appears 
to blast material away, or at least change its temperature enough to put the 
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Figure 6. Zoomed in view of single large oscillation region found in Figures [1] and [5] Panels 
(a) and (b) are from the 171A data, and panels (c) and (d) are derived from the 195A data (1 
July 1998). 

material out of the TRACE 171A passband. The physical phenomena captured 
by these observations imply that the resultant time-series have a significant back- 
ground trend which varies from location to location. The detrending timescale 
(500 seconds) removes the secular background intensity variations such as flaring 
and dimming. 

Results for the analysis in a three (120-240 seconds) and five minute frequency 
band (240-360 seconds) are given in Figure [HI and Figure [TO] rcpestively, and 
Table [3] Only three three-minute frequency band regions are found. Two of 
them (regions 1 and 3) are distant from the central active region over relatively 
dark pieces of corona where the signal is weak. Region 2 (the largest of the three) 
overlies a portion of the active region where many oscillations are detected in the 
five-minute frequency band; the detection in the three-minute frequency band 
may be due to multiple loops each oscillating in the five-minute frequency band 
that coincidcntally gives the appearance of a three-minute oscillation. 

The detection algorithm finds eight of the nine transversely oscillating as 
described in Aschwan den et aZ.I (|1999|) (their Figure 1, oscillations 1 through 9 
excepting 5). Many of the other claimed detections of Figure [TO] are in small 
regions distant from the flare site. These locations also typically show an oscilla- 
tion which is more of a gentle sway, that is, the oscillation decays very quickly. It 
is instructive to consider Figure EHe) in comparison to Figure EJb). The spatial 
probability distribution is much denser at all lengthscalcs in Figure EJc) than 
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Figure 7. Example image data and results from the analysis algorithm for TRACE 171 A 
data taken on 1 July 1998 (panel (a)). Panel (b) is the probability that a given pixel supports 
an oscillation within the frequency band indicated in panel (c), in this case, a five- minute 
oscillation frequency band. Also indicated in panels (a-c) are the detected oscillation regions. 
Panels (d — f) show similar maps to those shown and described in Figure [4] Table [2fc) shows 
the values for the detected region for the parameters listed in Table [T] 



in Figure [3jb) . Figure [TT] shows that the material of Figure [3je) has a greater 
proportion of higher-probability oscillations than the material in Figure E^b) . 
Given the differing spatial distribution, it is clear that everywhere in TRACE 
field of view on 14 July 1998 was much more likely to oscillate in the five-minute 
frequency band. Since the probability map is filtered for regions that show a 
locally high probability, many more regions are found in this dataset than in the 
previous dataset. 
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Figure 8. Example image data and results from the analysis algorithm for TRACE 195A 
data taken on 1 July 1998 (panel (a)) in the five-minute frequency band. Panels (b-f) show 
similar maps to those shown and described in Figure [7] Table [2{d) shows the values for the 
detected region for the parameters listed in Table [T1 



5. Detecting Other Oscillatory Signals with this Algorithm 

This paper is primarily concerned with detecting areas in the solar atmosphere 
that oscillate with a single frequency, described by Equation ([2}. Other oscilla- 
tions have been found in the solar atmosphere that are not perfectly described 
by Equation @. Brctthorst (1988, ch. 6.1.4) notes that Equation (0) still gives 
strongly peaked probability distributions close to the true frequencies even when 
the actual signal exhibits features not present in the model, such as periodic 
but non-harmonic oscillations and non-stationary and non-Gaussian noise. In 
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Figure 9. Example image data and results from the analysis algorithm for TRACE 171A 
data taken on 14 July 1998 (panel (a)) in the three-minute frequency band. Panels (b — f) 
show similar maps to those shown and described in Figure [7] Panels (b — f) show similar maps 
to those shown and described in [7] Table [3] shows the values for the detected region for the 
parameters listed in Table [T] 



the sections below we discuss some commonly occuring periodic signals in solar 
atmosphere, that contain features not modeled by Equation ([2]), and their effect 
on the probability distributions that are at the center of the proposed detection 
algorithm. 
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Figure 10. Example image data and results from the analysis algorithm for TRACE 171A 
data taken on 14 July 1998 (panel (a)) in the five-minute frequency band. Panels (b — f) show 
similar maps to those shown and described in Figure [7] Panels (b-f) show similar maps to 
those shown and described in Figure [7] Table [3] shows the values for the detected region for 
the parameters listed in Table [T] 



5.1. Decaying Oscillations 

Figure [T^] applies the single-frequency model to an example dataset based on 
the transverse loop oscillation described by Nakaria kov et aZ.I(|1999p . The data 
(shown in Figure EH a)) exhibits approximately the same properties as the true 
observation. Figure IT27 b) shows the measured frequency to be close to the true 
frequency of 3.9 mHz. For comparison, Figure [T^c) shows the same time-series 
as Figure 112( a) except with no decay, along with the probability distribution 
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Table 3. Quantities reported (see Tabled for the automatically detected regions found in 
14 July 1998 TRACE 171 A - see Figure [TO] for maps and plots of the detected regions. The 
quantities Q and E are also listed with their rank r when compared to all other regions 
found in the same dataset. Region 4 is the same region as that studied by Nakariakov et al. 



(1999) and Ireland and Dc Moortel (2002), and region 8 is the base of the same loop. The 
numbers Ax in the region column refer to the oscillations found manually in Figure 1 of 
lAschwanden et al. (1999). 
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for the frequency in Figure [T2Td) . The only difference between the two is the 
error in the determination of the frequency. This is because the decay in the first 
time-series decreases the signal to noise ratio, and so later portions of the time- 
series contribute less information to the determination of the frequency, and so 
the error increases. This shows that for decaying oscillations of the type already 
observed, the probability based methods of Section 2 yield good results. Further, 
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Figure 11. Distribution of non-zero probabilities for the 1 July and 14 July 17lA data. 
The dashed line refers to the 1 July 171A data analyzed in the three-minute frequency band 
(120-240 seconds), as shown in Figure (3jb). The solid line refers to the 14 July 17lA data 
analyzed in the five- minute frequency band (240-360 seconds), as shown in Figure 02e). 
Average values to each frequency distribution are indicated by the vertical lines of the same 
line style. See Section 14.21 for more detail. 



it also demonstrates that even although we used an inappropriate model, the 
model of Section 2 still gives a sufficiently peaked distribution at the right 
frequency to enable detection. 

5.2. Non-Stationary Frequency Oscillations 

Figure fTS] applies the single-frequency model to an example dataset based on 
the "tadpole" signature of Nakari akov et al] (|2004|) . The example data-set re- 
produces the extent of the signal as mentioned in INakariakov et aL\ (|2004|> . 
who claim that there is a time- varying (i.e., non-stationary) frequency present 
in a signal found by |Katsiyannis et al] (|2003p (see also I Williams et al\ 120011 
I Williams et al.[ I2002p in the Solar Eclipse Imaging System (SECIS) of Queens 
University, Northern Ireland eclispe data taken on 11 August 1999. Their obser- 
vation is modeled as having a signal to noise ratio of 1.5, and a non-stationary 
frequency that lowers by 5% over the range (Figure flBT a)). The result of an 
analysis using Equation (fT2"]) with Equation ^ (note that Equation (J3]), as 
used in the automated detection algorithm, is essentially the same) is shown in 
Figure ITBTb) . Figure ITBTc) shows the same extent of signal except that now the 
oscillation is stationary. Note that the algorithm will not detect the oscillation as 
being non-stationary or as not extending across the observation range, as neither 
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Figure 12. Effect of an exponential decay in the observed signal on the detection of an oscil- 
lation. Panel (a) shows the observed noisy time series signal (dotted line), the true (non- noisy) 
signal (solid line) and the estimated signal (dashed line). Panel (b) shows the probability-den- 
sity function for the observed signal as defined by Equation (£3j) - Panels (c) and (d) are the 
same as (a) except the exponential decay is no longer present. The probability-density function 
in Panel (d) is much narrower than in Panel (b). The effect of the exponential decay is to make 
the location of the average frequency more uncertain. This is because the signal-to-noise ratio 
decreases due to the exponential decay and therefore there is less information to determine the 
frequency when compared to the observed signal of panel (c). See Section l5.1l for more detail. 



of those effects are included in the model oscillation Equation @. However, 
the probability distribution functions Figures I13f b.d) are strongly and singly 
peaked, and they would therefore be detected given a wide enough detection 
window wi , u>2 (see Equation [7} ■ 

5.3. Multiple discrete frequencies 

The results stated in this paper are derived with regard to a single frequency 
model, Equation ([2]). Initial identification of a pixel as containing a significant 
oscillation relies on the integrated probability over a user-defined frequency band 
exceeding a high user-defined limit (in this case, 0.95, see Sectior J3~3]) . If there are 
multiple frequencies inside the user-defined frequency band and the integrated 
probability exceeds 0.95, then the algorithm will report a detection at that pixel 
and report the average frequency (and error) within that frequency band. It 
cannot do any more, since the model assumes the presence of only a single 
frequency. 
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Figure 13. Effect of a non-stationary frequency in the observed signal on the detection of 
an oscillation. Panel (a) shows the observed noisy time series signal (dotted line), the true 
(non- noisy) signal (solid line) and the estimated signal (dashed line). The true signal consists 
of noise at < t < 14 and a time-varying frequency at t > 14. Panel (b) shows the probabili- 
ty-density function for the observed signal as defined by Equation JSJ. Panels (c) and (d) show 
the same quantities as (a) and (b) respectively except that the frequency in the observed and 
true signals are constant at t > 14. The peak of the probability-density in panel (b) is at a 
lower frequency than in panel (d) because the observed signal contains more low frequencies. 
Note that in both cases the probability-density function is strongly peaked and so most of the 
probability lies in a narrow band around the peak. See Section l5.2l for more detail. 



Consider now the case of two well-separated frequencies such that one fre- 
quency lies inside the user-defined frequency band (wi, u>2, Equation Q) and the 
other lies outside the user-defined frequency band f Figure [Mk). Both frequencies 
are present in the observation, and both have the same amplitude and signal-to- 
noise ratio. The model however, supposes the presence of one frequency. Hence 
the probability will be split between these two frequencies (Figure ITlTb)). Which 
frequency has the highest probably depends on the evidence for each in the time- 
series, that is, the number of oscillations in the time series and the signal-to-noise 
ratio for each frequency present. In the case above, the most probable frequency 
is in the three-minute wave band. If one were looking in this frequency band, 
then a peak would be detected and the algorithm would decide if there were 
enough probability in the user-defined frequency band [wi, UJ2\ to claim a detec- 
tion. However if the user-defined frequency band [wi,W2] was in the five- minute 
frequency band, then since the majority of the probability lies outside this range, 
no detection would be claimed. However, we note that observations of multiple 
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Figure 14. Analysis of a time-series that has two widely separated present. In panel (a), the 
observed signal (dotted line) is a noisy observation of the true signal (solid). The estimated 
signal (dashed) is found using Equation 10, a single frequency model. Panel (b) shows the 
probability-density function (Equation (O) of the observed signal. In this case, the maximum 
probability density is close to the frequency of the higher frequency signal because the signal 
contains more information about this frequency than the lower frequency. See Section 15.31 for 
more detail. 



oscillations in single structures in the corona show one dominant oscillation and 
a second, much weaker oscillation (IVerwichte et al.l(|2004j) ). Therefore, the situa- 
tion of Figure Eta) - two oscillations of equal amplitude but different frequency 
- is observationally unlikely in the corona. However, one may use the approach 
here to first identify pixels that have a dominant oscillation present, and use other 
methods to examine for the presence of secondary, smaller amplitude oscillations. 
Deciding how many frequencies a time-series supports is specifically excluded 
by the model choice at the start, that is, a single oscillation plus noise. To prop- 
erly decide how many frequencies are present in a time-series, one must define a 
model for each case to be considered (cither no frequencies, one frequency, two 
frequencies, three frequencies, etcetera), calculate the probability of each model, 
and decide which model is the most probable. Bayesian probability tests hypothe- 
ses (given the data), and those hypotheses must be explicitly stated. Introducing 
other hypotheses, such as multiple frequency models, is asking a different ques- 
tion of the data. For example, multiple frequencies have long been observed 
in the type of sunspot observations examined by |Marsh, Ireland, and Kucera| 
(2008]), and therefore a single-oscillation model is inappropriate to start with; one 
simply would not use this model in this case. In addition, the approach taken by 
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|Marsh, Ireland, and Kucera| (|2008p is prohibitively expensive computationally 
at present, and so is not suitable for the automated detection of oscillations in 
large datasets. 



6. Conclusion 

Oscillations in the solar atmosphere have already demonstrated their worth as 
probes of the physical conditions present. SDO will be a major contributor to the 
study of solar atmospheric oscillations, and automated detection algorithms will 
be necessary in order to maximize the scientific potential here. The algorithm 
parameters in the filtering procedure (see Section [575]) will have to be tailored 
to SDO data; other algorithms will also probably have to undergo a "tweaking" 
phase to operate optimally. 

The algorithm described here is a first attempt at implementing an automated 
coronal oscillation detector based on a Bayesian understanding of probability. 
It shows promise in being able to find areas of the solar atmosphere that are 
highly likely to support an oscillation. Section @] shows that it is able to find both 
longitudinal and transversely oscillating loops at low estimated signal-to-noisc 
ratio over a complex background scene (the non-oscillating or slowly varying 
background, where slowly varying is understood as evolution longer than the 
characteristic timescale of the background trend subtraction). We note that even 
by eye, these events appear to be at low signal-to-noise ratio, and the estimates 
calculated here (Figures HUTU]) agree with this observation. 

We also examined the same quiet-Sun TRACE 1600 and I700A data as 
Mcintos h, Fleck, and Tarbell| (I2004p . In that paper, the authors calculate the 
travel time of waves in two different UV wave bands, with the understanding 
that the entire FOV contains oscillations. We find that almost the entire FOV 
contains oscillations. The algorithm returns positive detections over almost all of 
the FOV for other TRACE UV data (jMcAteer et all I20MJ) . In both cases the 
detection itself takes only a few seconds: however, the other parameters such as 
signal-to-noise ratio takes many multiples of the detection time. This suggests 
that a line production version which generates data products similar to Figures 
SI [TUJ or [5] be restricted to data where we know that only a small percentage of 
the FOV contains a signal, such as the higher temperature corona. An alternate 
mode of operation is that the reported quantities are restricted to detections 
only, with further quantities (such as signal to noise ratios) calculated offline by 
the interested user. 

The model oscillation that Equation @ describes, for a single pixel, a sinu- 
soidol intensity oscillation existing for the entire duration of the time-series, in 
the presence of Gaussian noise. Although this seems quite a restrictive model, it is 
clear from the results of Section 14721 that it is sufficient to enable the detection of 
decaying transverse oscillations. Indeed, llreland and De Moortell (|2002[) showed 
that the oscillation studied by [Nakariak ov et aL\ (j!999[) (region 4 in Table [3J 
displays non-Gaussian noise and a linear change in period over the measurable 
duration of the oscillatory motion, complications not contained in the model 
choice, Equation ^ . A transverse oscillation on a single pixel appears as a peri- 
odic non-sinusoidal change in intensity as the loop swipes across the field of view. 
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Despite the great difference between the model choice and the appearance of such 
an oscillation in this analysis, the resemblance between the two is sufficient for 
our algorithm to be able to detect transversely oscillating material. This is an 
example of the sufficiency of the model, (Equation @), to adequately describe 
a wide range of periodic behavior, as previously noted by [Bretthorst (1988). 

The algorithm runs quickly on the data (less than a minute, reported in Table 
2|) on an Apple MacBook (dual core, 1.6 Ghz processor), without any special 
efforts at algorithmic optimization. Analysis of binned (prepared) full-disk SDO 
data (see Table S]) is possible with existing computers. Note, however, that the 
analysis times quoted do not take into account the time required to prepare the 
data for automated oscillation analysis. This overhead will presumably be the 
same for all detection algorithms, and must be considered in design of a full-scale 
operational SDO oscillation-detection algorithm in order for it to run faster than 
the time taken to acquire the data. 

This paper introduces two new features to the discussion of automated detec- 
tion of oscillations in the solar atmosphere. Firstly, the algorithm is implemented 
using a Bayesian interpretation of probability as a "degree of belief" as opposed 
to the standard interpretation as "frequency of occurence" (leading to powerful 
and convenient formulae such as Equations [T3] and [3]) . At the core of the algo- 
rithm lies the probability that solar atmospheric time-series can be described as 
a single sinusoidal oscillation at a fixed frequency, subject to distortion Gaus- 
sian noise. As IBretthorstl (|1988p notes, this is a good approximation for many 
purposes. This algorithm does not say anything about the presence of two or 
more sinusoidal signals in a single time-series. However, it is certainly possible 
to develop an analysis algorithm to assign probabilities to each of the three 
hypotheses that the time time-series is either noise, contains a single frequency or 
contains two frequencies. Most of the required mathematics is quoted or referred 
to in this work; such an algorithm will be the described in future developments. 

Secondly, we have introduced "quality measures" in an attempt to grade the 
regions that survive the region finding and filtering process. This appears to 
be necessary given the large number of oscillating regions that can occur in 
a given datasct, for example the 14 July 1998 TRACE 17lA data. Further 
criteria may be set by individual users. For example high values of the ratio 
A0.95/A indicate that the region has a high proportion of pixels very probably 
supporting a frequency; a threshold could be set to filter only those regions that 
have a high proportion of very certain pixels. In addition, the current VO Event 
standard makes for the provision of extra algorithm information such as arbitrary 
parameters (R, m) to be carried along with any results. This means the user will 
be fully informed of all the parameter values used to obtain the results. 

Future algorithm development will concentrate on improving the probability 
map filtering (step 6 of Figure [2j and also Section 13. 3|) , extending the analysis 
to assign probabilities to multi-frequency models, and to distinguishing between 
longitudinal and transverse oscillations (it may be possible to assign a prob- 
ability that a given wave-mode has been observed). Understanding the wave 
mode demands an understanding of the structure on which it is supported, 
which naturally leads to the automated detection and characterization of loop 
structures, which is a complex topic by itself (see lA"schwand cn et aL\ (|2008p for 
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Table 4. Performance of the algorithm on the analyzed data compared to data 
acquisition time. Also shown is the projected performance on SDO data. 

data data size algorithm time [observation duration] 

1 July 256 X 256 30 [6231] 

201 samples 
31 s cadence 
14 July 368 x 368 25 [3431] 

47 samples 
31 s cadence 
raw SDO 4096 x 4096 7641 [2000] 

200 samples 
10 s cadence 
2x2x2 rebinned 2048 X 2048 955 [2000] 

SDO 200 samples 

20 s cadence 



a review). We hope that the algorithm presented here is a first step towards 
automated coronal seismology. 

Acknowledgements This work was supported by the NASA SESDA contract and a NASA 
Hcliophysics Guest Investigator award (NNG08EL33C). SOHO is a joint project of interna- 
tional co-operation by ESA and NASA. 



References 

Aschwanden, M.J., Fletcher. L., Schrijver, C.J., Alexander. D.: 1999, Coronal Loop Oscillations 

Observed with the Transition Region and Coronal Explorer. Astrophys. J. 520, 880-894. 

doi: 10. 1086/307502 
Aschwanden, M.J., de Pontieu, B., Schrijver, C.J., Title, A.M.: 2002, Transverse Oscillations 

in Coronal Loops Observed with TRACE II. Measurements of Geometric and Physical 

Parameters. Solar Phys. 206, 99-132. doi: 10. 1023/A: 1014916701283 
Aschwanden, M.J., Lee, J.K., Gary, G.A., Smith, M., Inhester, B.: 2008, Comparison of Five 

Numerical Codes for Automated Tracing of Coronal Loops. Solar Phys. 248, 359 — 377. 

doi: 10. 1007/sll207-007-9064-9 
Baycs, T.: 1763, An Essay towards solving a Problem in the Doctrine of Chances. Phil. Trans. 

Roy. Soc. 53, 370-418. http://rstl.royalsocietypublishing.org/content/53/370. 
Berger, T.E., de Pontieu, B., Fletcher, L., Schrij ver, C.J., Tarbell, T.D., Ti tle, A.M.: 1999, 

What is Moss? Solar Phys. 190, 409-418. doi 10.1023/A:1005286503963 
Bcrghmans, D., Clette, F.: 1999, Active region EUV transient brightcnings - First Results by 

EIT of SOHO JOP80. Solar Phys. 186, 207-229. 
Brctthorst, G.L.: 1988, Lecture Notes in Statistics: Bayesian Spectrum Analysis and Param- 
eter Estimation, Springer, Berlin Heidelberg. 
Chatficld, C: 1996, The analysis of time series: an introduction, 5th edn. Chapman and Hall, 

Florida. 
Cooley, J., Tukey, J.: 1965, An Algorithm for the Machine Computation of Complex Fourier 

Series . Math. Comp. 19, 297-301. 
De Moortel, I.: 2005, An overview of coronal seismology. Roy. Soc. London Phil. Trans. Series 

A 363, 2743-2760. 
Dc Moortel, I., Hood, A.W.: 2003, The damping of slow MHD waves in solar coronal magnetic 

fields. Astron. Astrophys. 408, 755-765. doi |10.1051/0004-6361:20030984| 



ts_preprint.tex; 7 July 2010; 0:36; p. 28 



28 Ireland et al. 

De Moortol, I., Hood, A.W.: 2004, The damping of slow MHD waves in solar coronal mag- 
netic fields. II. The effect of gravitational stratification and field line divergence. Astron. 

Astrophys. 415, 705-715. doi: 10.1051/0004-6361:20034233 
Dc Moortel, I., McAteer, R.T.J.: 2004, Waves and wavelets: An automated detection technique 

for solar oscillations. Solar Phys. 223, 1-11. doi |10.1007/sll 207-004-0806-7 
De Moortel, I., Ireland, J., Walsh, R.W.: 2000, Observation of oscillations in coronal loops. 

Astron. Astrophys. 355, L23-L26. 
Dc Moortel, I., Ireland, J., Hood, A.W., Walsh, R.W.: 2002, The detection of 3 and 

5 min period oscillations in coronal loops. Astron. Astrophys. 387, L13 — L16. 

doi: 10.1051/0004-6361:20020436 
de Pontieu, B., Erdclyi, R., De Moortel, I.: 2005, How to Channel Photosphcric Oscillations 

into the Corona. Astrophys. J. Lett. 624, L61-L64. doi 10.1086/430345] 
de Pontieu, B., Erdelyi, R., James, S.P.: 2004, Solar chromosphcric spicules from the leakage 

of photosphcric oscillations and flows. Nature 430, 536-539. doi 10.1038/nature02749 
de Pontieu, B., Berger, T.E., Schrijver, C.J., Title, A.M.: 1999, Dynamics of Transition Region 

'Moss' at high time resolution. Solar Phys. 190, 419-435. doi 10.1023/A:1005220606223 
Deforest, C.E., Gurman, J.B.: 1998, Observation of Quasi-periodic Compressive Waves in Solar 

Polar Plumes. Astrophys. J. Lett. 501, L217. doi 10.1086/311460 
Fletcher, L., de Pontieu, B.: 1999, Plasma Diagnostics of Transition Region "Moss" using 

SOHO/CDS and TRACE. Astrophys. J. Lett. 520, L135-L138. doi:10.108 6/312157| 
Grcchncv, V.V.: 2003, A method to analyze imaging radio data on solar flares. Solar Phys. 

213, 103-110. 
Gregory, P.C.: 2005, Bayesian Logical Data Analysis for the Physical Sciences: A Comparative 

Approach with 'Mathematica' Support, Cambridge University Press, Cambridge, UK. 
Ireland, J., Dc Moortel, I.: 2002, Application of wavelet analysis to transversal coronal loop 

oscillations. Astron. Astrophys. 391, 339-351. doi: 10.1051/0004-6361:20020643 
Ireland, J., Walsh, R.W., Harrison, R.A., Priest, E.R.: 1999, A wavelet analysis of active region 

oscillations. Astron. Astrophys. 347, 355-365. 
Jaynes, E.T.: 1987, Bayesian spectrum and chirp analysis. In: Smith, C.R., Erickson, G.J. (eds.) 

Maximum Entropy and Bayesian Spectral Analysis and Estimation Problems, 1 - 37. 
Jaynes, E.T.: 2003, Probability Theory: The Logic Of Science, 1st edn. Cambridge University 

Press, Cambridge, United Kingdom. 
Katsiyannis, A.C., Williams, D.R., McAteer, R.T.J., Gallagher, P.T., Keenan, F.P., Murtagh, 

F.: 2003, Eclipse observations of high-frequency oscillations in active region coronal loops. 

Astron. Astrophys. 406, 709-714. doi: 10.1051/0004-6361:20030458 
King, D.B., Nakariakov, V.M., Dcluca, E.E., Golub, L., McClements, K.G.: 2003, Propagating 

EUV disturbances in the Solar corona: Two-wavelength observations. Astron. Astrophys. 

404, L1-L4. doi: 10.1051/0004-6361:20030763 
Linncll-Nemec, A.F., Nemec, J.M.: 1985, A test of significance for periods derived using phasc- 

dispersion-minimization techniques. Astronom. J. 90, 2317-2320. doi: 10.1086/113936 
Marsh, M.S., Ireland, J., Kucera, T.: 2008, Bayesian Analysis of Solar Oscillations. Astrophys. 

J. 681, 672-679. doi: 10.1086/588751) 
McAteer, R.T.J., Gallagher, P.T., Bloomfield, D.S., Williams, D.R., Mathioudakis, M., 

Keenan, F.P. : 2004, Ultraviolet Oscillations in the Chromosphere of the Quiet Sun. 

Astrophys. J. 602, 436-445. doi: 10.1086/380835 
Mcintosh, S.W., de Pontieu, B., Tomczyk, S.: 2008, A Coherence-Based Approach for Tracking 

Waves in the Solar Corona. Solar Phys. 252, 321-348. doi:10.1007/sll207-008-9257-x 
Mcintosh, S.W., Fleck, B., Tarbell, T.D.: 2004, Chromosphcric Oscillations in an Equatorial 

Coronal Hole. Astrophys. J. Lett. 609, L95-L98. doi 10.1086/422748 
Nakari akov, V.M., King, D.B.: 2 007, Coronal Periodmaps. Solar Phys. 241, 397-409. 

doi: 10. 1007/sll207-007-0348-x 
Nakariakov, V.M., Vcrwichtc, E.: 2005, Coronal Waves and Oscillations. Living Reviews in 

Solar Physics 2, 3. 
Nakariakov, V.M., Ofman, L., Deluca, E.E., Roberts, B., Davila, J.M.: 1999, TRACE obser- 
vation of damped coronal loop oscillations: Implications for coronal heating. Science 285, 

862-864. doi: 10.1126/science. 285. 5429. 862 
Nakariakov, V.M., Arbor, T.D., Ault, C.E., Katsiyannis, A.C., Williams, D.R., Keenan, F.P.: 

2004, Time signatures of impulsively generated coronal fast wave trains. Mon. Not. Roy. 

Astron. Soc. 349, 705-709. doi |10.1111/j. 1365 -2966. 2004. 0753Txl 



ts_preprint.tex; 7 July 2010; 0:36; p. 29 



Automated Detection of Oscillating Regions in the Solar Atmosphere 29 

Nightingale, R.W., Aschwanden, M.J., Hurlburt, N.E.: 1999, Time Variability of EUV 

Brightcnings in Coronal Loops Observed with TRACE. Solar Phys. 190, 249-265. 

doi: 10. 1023/A: 1005211618498 
Ofman, L., Romoli, M., Poletto, C, Noci, G., Kohl, J.L.: 1997, Ultraviolet Coronagraph 

Spectrometer Observations of Density Fluctuations in the Solar Wind. Astrophys. J. 

Lett. 491, Llll. doi 10.1086/311067 
Ruanaidh, J.K., Fitzgerald, W.J.: 1996, Numerical Bayesian methods applied to signal 

processing, Springer, New York/London. 
O'Shea, E., Banerjee, D., Doyle, J.G., Fleck, B., Murtagh, F.: 2001, Active region oscillations. 

Astron. Astrophys. 368, 1095-1107. doi |10.105l7 0004-6361:20010073"1 
Robbrecht, E., Verwichte, E., Berghmans, D., Hochcdez, J.F., Poedts, S., Nakariakov, V.M.: 

2001, Slow magnctoacoustic waves in coronal loops: EIT and TRACE. Astron. Astrophys. 

370, 591-601. doi 10.1051/0004-6361:20010226 
Schuster, A.: 1898, On the Investigation of Hidden Periodicities with Application to a Supposed 

26 Day Period of Meteorological Phenomena. Terr. Mag. 3, 13-13. 
Sych, R.A., Nakariakov, V.M.: 2008, The Pixelised Wavelet Filtering Method to Study Waves 

and Oscillations in Time Sequences of Solar Atmospheric Images. Solar Phys. 248, 395- 

408. doi 10.1007/sll207-007-9005-7 
Torrencc, O, Compo, G.P.: 1998, A Practical Guide to Wavelet Analysis. Bull. Am. Meteor. 

Soc. 79, 61-78. doi. 10.1175/1520-0477(1998)079 
Verwichte, E., Nakariakov, V.M., Ofman, L., Deluca, E.E.: 2004, Characteristics 

of transverse oscillations in a coronal loop arcade. Solar Phys. 223, 77-94. 

doi:10.1007/sll207-004-0807-6 
Williams, D.R., Phillips, K.J.H., Rudawy, P., Mathioudakis, M., Gallagher, P.T., O'Shea, 

E., Keenan, F.P., Read, P., Rompolt, B.: 2001, High-frequency oscillations in a 

solar active region coronal loop. Mon. Not. Roy. Astron. Soc. 326, 428-436. 

doi: 10.1046/j. 1365-8711. 2001. 04491.x 
Williams, D.R., Mathioudakis, M., Gallagher, P.T., Phillips, K.J.H., McAteer, R.T.J. , Keenan, 

F.P., Rudawy, P., Katsiyannis, A.C.: 2002, An observational study of a magncto- 
acoustic wave in the solar co rona. Mon. Not. Roy. Astron. Soc. 336, 747-752. 

doi |l03046/j.l365-8711.2002.05764~x1 

Appendix 

A. The General Linear Model 

A time-series is a special case of a more g eneral linear model descripti on. In 



this section we briefly recap the argument of O Ruanaidh and Fitzgerald (I1996[) 
in deriving Equation §5§. In this description, a signal d(i) observed at times ij 
(1 < i < N), is modeled as 



M 

d(i) = ^2 b k9k(i) + x{i), (8) 

fc=i 

with M basis functions <?&(«), each of amplitude bk, evaluated at time ti (param- 
eterized by {uj}) and Gaussian distributed noise x(i) of mean and standard 
deviation a. In the context of Section [^ the observation D is equivalent to 
the signal d(i), and the hypothesis H is that the data can be described by the 
right hand side of Equation ((5J , including the noise. In matrix form, the above 
equation may be written as 

d = Gb + x, (9) 
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where d is a N x 1 matrix of the observed data, and x is a N x 1 matrix of 
identically distributed and independent Gaussian noise samples. The matrix G 
is size N x M; each column of G is one of the basis functions evaluated at all ij. 
The matrix b is a M X 1 matrix, the linear coefficients of each of the (column) 
basis functions in G. The likelihood function of the observed data is 



p(d\{u},<r,b,I) = 



1 


r t ' 

X X 

2a 2 _ 


Gb) T (d- 




1 


Gb)" 


(2™ 2 )^ 2 " P 




2a 2 





(10) 



where {w} parameterize the basis functions g^ and hence G. This equation is 
derived by multiplying together the probability distributions of the noise x% at 
each time i. The exponent shows that maximizing the probability is equivalent to 
minimizing the difference between the data and the basis functions (weighted by 
their amplitudes); indeed, this equation forms the basis of least-squares fitting 
in the presence of Gaussian noise. 

In analysis, one is primarily interested in the values of the parameters Jul, and 
secondarily interested in the other values such as b and a. O Ruanaidh and Fitzgerald] 
(|1996p and IBretthorstl (|1988p describe the process by which the "nuisance pa- 
rameters" b and a are removed from further consideration by marginalization. 
The statement of Bayes' theorem for the general linear model, using Equations 

© and dinD is 



p({u},b,<T\d,I) 



p({o;},b,<x|J)p(d|{a;},<7,b,J) 

p(d\I) 



Integrating over b and a using a prior p({uj}, b, cr\I), removes these variables 
from further explicit consideration, and is an example of Bayesian marginal- 
ization. On integration, this obtains the marginal posterior distribution for the 
parameters {uj}: 



*>({«}|d,/) 



/ [ p{{cj},b,<T\d,I)dadb. 



(11) 



6 Ruanaidh and Fitzgerald (|1996p and IBretthorst (1988) use uniform priors for 
the amplitude parameters b (p(b, I) = constant) and the Jeffreys prior (p(a\I) oc 
1/<t) for a. On integration, 



p({w}|,d,J)« 



d T d - d T G G T G 



G T d 



{M-N)/2 



y/Act (G T G) 



(12) 



Equation (TT2l is a function of {w} only; the standard deviation (i.e., the noise 
level of the time-series) or the amplitude of the basis functions (the values b) need 
not be known in order for estimates of {w} to be found. This is a very powerful 
equation, with clear application to solar time-series analysis, where estimates of 
the noise level and oscillation amplitude in the data are often difficult to obtain 
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by direct fitting of model functions to the observed time-series. It should also be 
noted that Equation (fT12j) arises from Equations ^ and (fTU|) . which is the basis 
of a general least squares fit to the data given Gaussian distributed noise. 



B. Estimating Basis Function Amplitudes and Variance 

Maximizing Equation (|10|) with respect to b (on substitution of Equation ((9])) 
leads to an amplitude estimate: 

b=(G T G) _1 G T d. (13) 

This maximization is identical to the "least squares" fit to the data for a given 
value of {uj}. Given these amplitudes, the model fit is then 

f = Gb. (14) 



In addition, O Ruanaidh and Fitzgerald <|1996|) show that an estimate (found by 

dizing the amplitudes) to the Gaussian 

[d T d - f T f] (15) 



maximizing the posterior after marginalizing the amplitudes) to the Gaussian 
variance is 

1 



a 2 



N-M 



This estimated variance is the data energy minus the estimated signal energy, 
divided by the number of degrees of freedom. Note that b and a 2 arc functions 
of the analyzing frequency [to] . Equations ([T3")) and (TT5T) arc calculated for the 
oscillating regions detected via Section [3] and returned as part of the results in 
Section E] (see Figures HI [5j El and EJ. 
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